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A model including two nonlinear chains with linear and nonlinear couplings between them, and 
opposite signs of the discrete diffraction inside the chains, is introduced. In the case of the cubic 
[x'^'] nonlinearity, the model finds two different interpretations in terms of optical waveguide arrays, 
based on the diffraction-management concept. A continuum limit of the model is tantamount to a 
dual-core nonlinear optical fiber with opposite signs of dispersions in the two cores. Simultaneously, 
the system is equivalent to a formal discretization of the standard model of nonlinear optical fibers 
equipped with the Bragg grating. A straightforward discrete second-harmonic-generation [x''^'] 
model, with opposite signs of the diffraction at the fundamental and second harmonics, is introduced 
too. Starting from the anti-continuum (AC) limit, soliton solutions in the x''^' model are found, 
both above the phonon band and inside the gap. Solitons above the gap may be stable as long as 
they exist, but in the transition to the continuum limit they inevitably disappear. On the contrary, 
solitons inside the gap persist all the way up to the continuum limit. In the zero-mismatch case, 
they lose their stability long before reaching the continuum limit, but finite mismatch can have a 
stabilizing effect on them. A special procedure is developed to find discrete counterparts of the 
Bragg-grating gap solitons. It is concluded that they exist all the values of the coupling constant, 
but are stable only in the AC and continuum limits. Solitons are also found in the x'^' model. They 
start as stable solutions, but then lose their stability. Direct numerical simulations in the cases of 
instability reveal a variety of scenarios, including spontaneous transformation of the solitons into 
breather- like states, destruction of one of the components (in favor of the other), and symmetry- 
breaking effects. Quasi-periodic, as well as more complex, time dependences of the soliton amplitudes 
are also observed as a result of the instability development. 



I. INTRODUCTION 

A. Objectives of the work 

Solitary-wave excitations in discrete nonlinear dynam- 
ical models (lattices) is a subject of great current inter- 
est, which was strongly bolstered by experimental ob- 
servation of solitons in arrays of linearly coupled opti- 
cal waveguides [1] and development of the diffraction 
management (DM) technique, which makes it possible 
to effectively control the discrete diffraction in the ar- 
ray, including a possibility to reverse its sign (make the 
diffraction anomalous) [2,3]. It has recently been shown 
that a lattice subject to periodically modulated DM can 
also support stable solitons, both single-component ones 
[4,5] and two-component solitons with nonlinear coupling 
between the components via the cross-phase-modulation 
(XPM) [6]. 

Two-component nonlinear-wave systems, both contin- 
uum and discrete, which feature a linear coupling be- 
tween the components, constitute a class of media which 
can support gap solitons (GSs). A commonly known ex- 
ample of a continuum medium that gives rise to GSs is 
a nonhnear optical fiber carrying a Bragg grating [7,8], 
whose standard model is based on the equations 
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where 5'(x,t) and ^{x,t) are amplitudes of the right- 
and left-propagating waves, and the Bragg-reflection co- 
efficient is normalized to be I. Another optical system 
that may give rise to GSs is a dual-core optical fiber 
with asymmetric cores, in which the dispersion coeffi- 
cients have opposite signs [9]. 

In this work, we demonstrate that the use of the DM 
technique provides for an opportunity to build a dou- 
ble lattice in which two discrete subsystems with oppo- 
site signs of the effective diffraction are linearly coupled, 
thus opening a way to theoretical and experimental study 
of discrete GSs, as well as of solitons of different types 
(solitons in linearly coupled lattices with identical dis- 
crete diffraction in the two subsystems have recently been 
considered in Ref. [10]; a possibility of the existence of 
discrete GSs in a model of a nonlinear-waveguide array 
consisting of alternating cores with two different values 
of the propagation constant was also considered recently 
[II]). The objective of the work is to introduce this class 
of systems and find fundamental solitons in them, in- 
cluding the investigation of their stability. We will also 
consider, in a more concise form, another physically rele- 
vant possibility, viz., a discrete system with a second- 
harmonic-generating (SHG) nonlinearity, in which the 
diffraction has opposite signs at the fundamental and sec- 
ond harmonics. Solitons will be found and investigated 
in the latter system too. 
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It is relevant to start with equations on which our x^^^ 
model (the one with the cubic nonlinearity) is based, 

i^ = - (C^2 + ?) Vn - {\^nf + f3 Icl^nf) i>n ' Kcf>n = 0, 
= S ■ (CA2 +q)cj)n- (|<^n|' + f3 l^nl') <^n " ^V'n = 0, 

where tpnit) and ^„(t) are complex dynamical variables 
in the two arrays (sublattices), k and /3 being coeffi- 
cients of the linear and XPM coupling between them, 
and t is actually not time, but the propagation distance 
along waveguides, in the case of the most physically rel- 
evant optical interpretation of the model. The operators 
CA2-0„ = C(7/;„+i +-0,1-1 - 2V'„) and {-6)CA2<j)n = 
{—6) C {4>n+i + 4'n-i ^ '^4>n) represent discrete diffrac- 
tion induced by the linear coupling between waveguides 
inside each array, the diffraction being normal in the first 
sublattice and anomalous in the second, with a negative 
relative diffraction coefficient —5 and intersite coupling 
constant C (one may always set C > 0, which we as- 
sume below). Physical reasons for having — 5 < are 
explained below. Finally, the real coefficient q accounts 
for a wavenumber mismatch between the sublattices. 

We also choose a similar SHG model, following the 
well-known pattern of discrete SHG systems with nor- 
mal diffraction at both harmonics [12], [13]: 

= -CA2V„ - CAn, (4) 
= 5CA20„ - 02 _ (5) 

where the asterisk stands for the complex conjugation 
and K is a real mismatch parameter. In this case too, we 
assume —5<Q. 

There are at least two different physical realizations of 
the x^^^ model based on Eqs. (2) and (3). First, one may 
consider two parallel arrays of nonlinear waveguides with 
different effective values n!^^^ and rS'^^ of the refractive in- 
dex in them corresponding to a given (oblique) direction 
of the light propagation. To this end, the waveguides 
belonging to the two arrays may be fabricated from dif- 
ferent materials; alternatively, they may simply differ by 
the transverse size of waveguiding cores, or by the refrac- 
tive index of the filling between the cores, sec e.g.. Fig. 1. 
The difference in the effective refractive index gives rise 
to the mismatch parameter q in Eqs. (2) and (3). More 
importantly, it may also give rise to different coefficients 
of the discrete diffraction. Indeed, the DM technique 
assumes launching light into the array obliquely, the ef- 
fective diffraction coefficient in each array being [2] 

£>(i'2) =2(7^2 cos (fc^'-'^rf), (6) 

where d is the spacing of both arrays, and A;^'^^ are trans- 
verse components of the two optical wave vectors. As it 



follows from Eq. (6), the diffraction coefficients are dif- 
ferent li k^l'> ^ kf'. 

Despite the fact that fc^^-* and k^f' are assumed dif- 
^cknt, we assume that the propagation directions of the 
light beams are parallel in the two arrays, as a conspicu- 
(^&^ walkoff (misalignment) between them will easily de- 
stroy any coherent pattern. On the other hand, the light 
coupled into both arrays has the same frequency, hence 
the absolute values of the two wave vectors are related 
as follows: k^^^ /k'^'^^ = n''^yn^'^\ where n^-'^'^' are the 
above-mentioned effective refractive indices. Combining 
the latter relation and the classical refraction law, and 
taking into regard the condition that the propagation di- 
rections arc parallel inside the arrays, one readily arrives 
at the conclusion that 

fc«/fcf =n«/n(2). (7) 

Note that the two incidence angles 9^^''^'> (at the inter- 
face between the arrays and air) are related in a similar 
way, (sin^^^^) / = n'^^^/n'^^\ hence the incident 

beams (in air) must be misaligned, in order to be aligned 
in the arrays. 




FIG. 1. Two parallel asymmetric arrays of optical waveg- 
uides that are described by Eqs. (2) and (3), provided that 
parallel beams propagate obliquely across both arrays. Ar- 
rows indicate misaligned directions at which light is coupled 
into the arrays; inside them, both propagation directions are 
identical. 

Equation (6) shows that there is a critical direction of 

(12,) 

the beam in each array, corresponding to k]^' 'd — tt/2, 
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at which the effective diffraction cocfhcicnt changes its 
sign [2]. Due to the difference between fc^' and , the 
critical directions are different in the two arrays. Then, if 
the common propagation direction in the arrays is chosen 
to be between the two critical directions, Eq. (6) gives 
different signs of the two diffractive coefficients. Note 
that this interpretation of the model implies no XPM 
coupling between the arrays, i.e., /3 = in Eqs. (2) and 
(3). 

An alternative realization is possible in a single array of 
himodal optical fibers, into which two parallel beams with 
orthogonal polarizations, u and v, are launched obliquely. 
If the two polarizations are circular ones, then /3 = 2 in 
Eqs. (2) and (3), and the asymmetry between the beams, 
which makes it possible to have different signs of the co- 
efficient (6) for them, may be induced by birefringence, 
which, in turn, can be easily generated by twist applied 
to the fibers [14] . The birefringence also gives rise to the 
mismatch q. As for the linear mixing between the two 
polarizations, which is assumed in the model, it can be 
easily induced if the fibers are, additionally, slightly de- 
formed, having an elliptic cross section [14]. If the two 
polarizations are linear, then the birefringence is induced 
by the elliptic deformation, and the linear mixing is in- 
duced by the twist, the XPM coefficient being 2/3 in this 
case (assuming that, as usual, the birefringence makes 
it possible to neglect four-wave mixing nonlinear terms 
[14]). 

It is interesting to note that the discrete model based 
on Eqs. (2) and (3) with k = is exactly tantamount 
to a formal discretization of the above-mentioned contin- 
uum model which was introduced in Ref. [9] to describe 
a dual-core optical fiber with opposite signs of disper- 
sion in the cores. Another quite noteworthy feature of 
the present model is that, if /3 = 2, it turns out to be 
formally equivalent to a discretization of the standard 
Bragg-grating model (1), which is produced by replacing 
^ - *„_i) /2 and ^ - $„_i) /2. 

Indeed, making the substitution ("staggering transfor- 
mation" ) 

*„ = i>„,$„ = i"Vn, (8) 

one concludes that the discrete version of Eqs. (1) takes 
precisely the form of Eqs. (2) and (3) with 5 = 1, q = 
2C, K = 1, and /3 = 2. 

B. The linear spectrum 



inside the gap of the system's linear spectrum one may 
find discrete gap solitons, i.e., counterparts of the GSs 
found in the continuum version of the model in Ref. [9] . 
Secondly, solitons specific to the discrete model may be 
found above the phonon band. To analyze these possibil- 
ities, an asymptotic expression for the tail, 

V'n,^n ~ exp(iAt- A|n|) (9) 

is to be substituted into the linearized version of Eqs. (2) 
and (3). 

Investigating the possibility of the existence of solitons 
above the phonon gap, it is sufficient to focus on the 
particular case 6=1 and q = 0, when the system's spec- 
trum takes a simple form (we have also considered more 
general cases with positive 5 different from 1 and q ^ 0, 
concluding that they do not yield anything essentially 
different from this case). The final result, produced by 
a straightforward algebra, is that solitons are possible in 
the region 

A^ > A^dg, = 16C^ + k\ (10) 



ibAedge being edges of the phonon band. In what follows 
below, we will assume A > 0, as in this case positive and 
negative values of A are equivalent. 

To understand the possibility of the existence of the 
discrete GSs, we, first, set (5 = 1 as above, but keep the 
mismatch q as an arbitrary parameter. Then, the gap is 
easily found to be 

r q^ + K'^ if qjO, 

k^<Kl,^=\ ifO<(7<4C, (11) 

[ (g - ACf + if q^4C 

(recall that, by definition, C > 0). An essential role of 
the mismatch parameter is that it makes the gap broader 
if it is negative. 

In the more general case, 5^1, two different layers 
can be identified in the gap, similar to what was found in 
the continuum limit [9]. For instance, if g = 0, the inner 
and outer layers are 

< A^ < — ^^K^, and — ^^k^ < A^ < (12) 
(,5+1)' {6+lf 



Before proceeding to the presentation of numerical re- 
sults for solitons found in the system of Eqs. (2) and (3), 
it is relevant to understand at which values of the propa- 
gation constant A (spatial frequency) solitons with expo- 
nentially decaying tails may exist in this model. There 
are two regions in which they may be found. Firstly, 



(in the case 5=1, the outer layer disappears). The 
difference between the layers is the same as in the con- 
tinuum limit [9] : in the outer layer, solitons, if any, have 
monotonically decaying tails, i.e., real A in Eq. (9), while 
in the inner layer A is complex, and, accordingly, soliton 
tails are expected to decay with oscillations. 
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C. The structure of the work 

The rest of the paper is organized as follows. In sec- 
tion II, we display results for solitons found above the 
phonon band, i.e., in the region (10). The evolution of the 
solitons is monitored, starting from the anti-continuum 
(AC) limit C = 0, and gradually increasing C. Any 
branch of soliton solutions in this region must disappear, 
approaching the continuum limit. Indeed, as the radia- 
tion band (frequently called "phonon band" , referring to 
linear phonon modes in the lattice dynamics) becomes 
infinitely broad in this limit, see Eq. (10), the solution 
branch with A = const will crash hitting the swelling 
phonon band. However, in many cases the soliton of this 
type is found to remain stable as long as it exists, so 
it may be easily observed experimentally in the optical 
array. 

In section III we present results for solitons existing 
inside the gap. In the outer layer [which is defined as per 
Eq. (12), provided that 5 7^ 1], we were able to find only 
solitons of an "antidark" type, that sat on top of a non- 
vanishing background. However, in the inner layer [recall 
it occupies the entire gap in the case 5=1, according to 
Eq. (12)], true solitons are easily found (in accord with 
the prediction, their tails decay with oscillations). In the 
case q = 0, these solutions appear as stable ones in the 
AC limit, get destabilized at some finite critical value of 
C, and continue, as unstable solutions, all the way up to 
the continuum limit, never disappearing. It is quite inter- 
esting that sufficiently large negative mismatch strongly 
extends the stability range for these solitons. 

As was mentioned above, the x*-^-* rnodel based on Eqs. 
(2) and (3) may be considered as a discretization of the 
standard gap-soliton system (1). In this connection, it is 
natural to search for discrete counterparts of the usual 
GSs in the latter system. However, the discrete GSs 
found in section III do not have any counterpart in the 
continuum system (1), as the staggering transformation 
(8) makes direct transition from the discrete equations 
(2) and (3) to the continuum system (1) impossible. At 
the end of section III, we specially consider discrete soli- 
tons which are directly related to GSs in the system (1). 
We find that such solitons exist indeed at all the values 
of C, their drastic difference from those found in sections 
II and III is that they are essentially complex solutions 
to the stationary version of Eqs. (2) and (3). At all fi- 
nite values of C, they are unstable, but the instability 
asymptotically vanishes in the AC and continuum limits, 
C ^ and C ^ 00. 

In section IV, we briefly consider the SHG model (4), 
(5). Solitons are found in this model too, and their sta- 
bility is investigated. When the solitons are linearly un- 
stable, the development of their instability is examined 
(in all the sections II, III, and IV) by means of direct nu- 
merical simulations, which show that the instability may 



initiate a transition to a localized breather, or to lattice 
turbulence, or, sometimes, complete decay of the soliton 
into lattice phonon waves. 

II. SOLITONS ABOVE THE PHONON BAND 

A. General considerations 

Stationary solutions to Eqs. (2)- (3) are sought for the 
form 

V„ = e*^*w„, 0„ = e'^'vn , (13) 

where A is the propagation constant defined above. In 
figures displayed below, the stationary solutions will be 
characterized by the norms of their two components, 

+00 +00 

p^^ J2 ui J2 <■ (14) 

n— — oo n— — oo 

Once such solutions are numerically identified by means 
of a Newton-type numerical scheme, we then proceed to 
investigate their stability, assuming that the solution is 
perturbed as follows: 

V'n = [un + eOn exp(2a;t) -I- ebn exp(— ia;*i)] exp(iAt), (15) 
0n = [vn + ec„ exp(zciji) -I- dn exp(— icij*t)] cxp(iAt), (16) 

where e is an infinitesimal amplitude of the perturba- 
tion, and u) is the eigenvalue corresponding to the linear 
(in)stability mode. The set of the resulting linearized 
equations for the perturbations {a,b*,c,d*;u)} is subse- 
quently solved as an eigenvalue problem. This is done by 
using standard numerical linear algebra subroutines built 
into mathematical software packages [15] . If all the eigen- 
values u) are purely real, the solution is marginally stable; 
on the contrary, the presence of a nonzero imaginary part 
of uj indicates that the soliton is unstable. When the so- 
lutions were unstable, their dynamical evolution was fol- 
lowed by means of fourth-order Runge-Kutta numerical 
integrators, to identify the development and outcome of 
the corresponding instabilities. 

In what follows below, we describe different classes of 
soliton solutions, which are generated, in the AC limit, 
by expressions with diff'erent symmetries. Still another 
class of solitons, which carries over into the usual GSs in 
the continuum system (1), will be considered in the next 
section. 

B. Solution families which are symmetric in the 
anti-continuum limit 

As it was said above, in this section we set 5 = 1 and 
q = 0, since comparison with more general numerically 
found results has demonstrated that this case adequately 
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represents the general situation, as concerns the existence 
and stabiUty of soUtons. Figure 2 shows a family of soli- 
ton solutions found for k = 0.1, A = 2 and /? = 0, as 
a function of the coupling constant C . In this case, the 
family starts, in the AC hmit (C = 0), with a solution 
that consists of a symmetric excitation localized at a sin- 
gle lattice site no, with 



/A- 



1 + ' 



(17) 



and terminates at finite C. Figure 2 demonstrates that 
this branch is always unstable. The termination of the 
branch happens when it comes close to the phonon band, 
that swells with the increase of C. The branch termi- 
nates at C = 0.464, when the upper edge of the band is 
at Aedgc = Vk^ + 16C2 w 1.859, according to Eq. (10). 
This value is still smaller than the fixed value of the soli- 
ton's propagation constant, A = 2, for which the soliton 
branch is displayed in Fig. 2. The branch, if it could be 
continued, would crash into the upper edge of the phonon 
band at C = 0.499. The slightly premature termination 
of this soliton family is a consequence of the nonlinear 
character of the solutions, as the above prediction for 
the termination point was based on the linear approxi- 
mation. 





FIG. 2. The top panel shows the norms Pu (lower curve) 
and Pv (upper curve) of the two components of the soliton 
solution vs. C, up to the point where the branch terminates. 
The next set of panels shows two examples of the solution at 
C = 0.1 (the upper row) and C = 0.464 (just near the ter- 
mination point of the branch; the lower row), together with 
the spectral planes of the corresponding linear stability eigen- 
values (the vertical and horizontal coordinates in the plane 
correspond to the imaginary and real parts of u)). The pro- 
files of the Un and Vn components are shown, respectively, by 
circles and stars. These solutions are always unstable. The 
bottom panel shows the imaginary part of the single unstable 
eigenfrequency vs. C. 



An example of the development of the instability of 
this solution, as found from direct simulations of the full 
equations (2) and (3), is given in Fig. 3 for C = 0.1. It is 
clearly seen that the unstable soliton turns into a stable 
breather. 
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FIG. 3. Evolution of the unstable soliton from Fig. 2) in 
the case C = 0.1, /3 = 0, and k = 0.1. The top panel shows the 
fields' spatial profiles (the circles correspond to ItAnP, and the 
stars to If/inp) for f = 4 (left panel), t — 192 (middle panel) 
and t = 196 (right panel). The first profile is nearly identi- 
cal to the initial condition, while the other two were chosen 
close to points where the oscillating amplitude of the resul- 
tant breather attains its maximum and minimum. The bot- 
tom panel shows the field evolution at the central lattice site 
{n = 50), clearly demonstrating the breathing nature of the 
established state. The solid and dashed lines are, respecively, 
I ■(/'sop and |05op. In this case, the instability growth rate of 
the initial soliton is ~ 0.8; in view of this large value, it was 
not necessary to add any initial perturbation to trigger the 
instability. 



C. Solution families which are anti-symmetric in the 
anti-continuum limit 



Another branch of solutions is initiated, in the AC 
limit, by an anti-symmetric excitation localized at a sin- 
gle lattice site, cf. Eq. (18): 



/A- 



On the contrary, in the presence of XPM with the phys- 
ically relevant value of /3 = 2, a similar solution branch, 
found for the same values k = 0.1 and A = 2, is stable 
for all C, until it terminates at C = 0.499. Note that, 
at this point, the upper edge (10) of the phonon band is 
Aedge = 1.999. which is extremely close to A = 2, i.e., the 
termination of the solution family is indeed accounted for 
by its crash into the swelling phonon band. Details of this 
stable branch are shown in Fig. 4. 



Direct simulations of this solution have corroborated 
its stability (details are not shown here). In fact, in all 
the cases when solitons are found to be stable in terms 
of the linearization eigenvalues (sec other cases below), 
direct simulations fully confirm their dynamical stability. 




FIG. 4. The same as in Fig. 2, but for the case (3 = 2. 
The middle and lower panels display examples of the soliton 
solutions at C = 0.1 and C = 0.499, respectively. In the top 
panel, the upper and lower curves now correspond to the Vn 
and Un components, i.e., opposite to the case shown in Fig. 
1. Notice that this branch is always stable until it termi- 
nates, therefore the figure does not contain a counterpart of 
the dependence shown in the bottom panel of Fig. 2. 



1 + /? 



(18) 



The solution belonging to this branch is shown in Fig. 
5 for the same values of parameters as in Fig. 2, i.e., 
5 = 1, K ~ 0.1, A = 2, and (3 = 0. With the increase 
of C, this branch picks up an oscillatory instability at 
C w 0.257, and terminates at (7 « 0.407. UnUke the 
solutions displayed above, the termination of this branch 
occurs not through its crash into the phonon band, but 
via a saddle-node bifurcation. The latter bifurcation im- 
plies a collision with another branch of solutions. That 
additional branch (which is strongly unstable) was found 
but is not shown here. 

In fact, the numerical algorithm is able to capture 
other solutions (unstable ones) past the point C « 0.407 
at which the present solution terminates. The newly 
found solutions are shown in the bottom part of Fig. 
5. However, the new family cannot be continued beyond 
C = 0.467 [cf. the termination point C = 0.464 for the 
solutions initiated in the AC limit by the expression (17)]. 
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FIG. 5. The top panel shows the branch of solutions start- 
ing from the expression (18) in the AC limit. Two particular 
examples arc shown for C = 0.1 (a stable solution, the upper 
row) and C = 0.257 (at the onset of the oscillatory instability, 
the lower row). The circles and stars again denote the u„ and 
v„ components, respectively. The panel beneath this displays 
the instability growth rate, Imw, vs. C. Finally, the bottom 
panels show the profiles and linear stability eigenvalues for 
another solution, found beyond the termination point of this 
branch at C = 0.407. Two examples of the new solution are 
given for C = 0.407 (the upper row) and C = 0.467 (the lower 
row). These two points are very close to the beginning and 
termination of the new branch). 

The development of the oscillatory instability of the 
solution shown in Fig. 5 was also studied in direct simu- 
lations. It leads to onset of a state where one component 
of the soliton is fully destroyed [it cannot completely dis- 
appear, due to the presence of the linear couplings in 
Eqs. (2) and (3), but it is reduced to a level of small 
random noise] . An example of this is given in Fig. 6 for 
C = 0.3, /3 = and K = 0.1. The instability (with the 
initial growth rate 0.07 in this case) develops after t « 60, 
destroying one component of the soliton in favor of fur- 
ther growth of the other one. In this uniformly 



distributed noise perturbation of an amplitude lO"** was 
added to accelerate the onset of the instability, as the 
initial instability is very weak (which implies that the 
unstable soliton may be observed in experiment). 




45 50 55 30 40 50 60 70 

n n 




t 

FIG. 6. Dynamical development of the oscillatory instabil- 
ity of the anti-symmetric solution for C = 0.3, /3 = and 
K = 0.1. The meaning of the symbols is as in Fig. 3. The top 
left and right panels show the field configuration at t = 4 and 
t = 196, respectively. The bottom panel once again shows the 
field evolution at the central site. 



A counterpart of the solution from Fig. 5, but with 
(3 — 2, rather than /3 = 0, is shown in Fig. 7. This 
branch is always unstable (i.e., in the case of the solu- 
tions starting from the anti-symmetric expression in the 
AC limit, the XPM nonlinearity destabilizes the solitons, 
while in the case of the branch that was initiated by the 
symmetric expression in the AC limit, the same XPM 
nonlinearity was stabilizing). It terminates at C « 0.219, 
again through a saddle-node bifurcation. As in the pre- 
vious case, a new family of solutions can be captured 
by the numerical algorithm past the termination point. 
The new family is found for 0.22 < C < 0.498, and it is 
also shown in Fig. 7. Comparing the value Aedge = 1-995 
given by Eq. (10) in this case with the actual value A = 2 
of the soliton's propagation constant, we conclude that 
the termination of the latter branch is caused by its col- 
lision with the phonon band. Notice also that the latter 
branch becomes unstable only very close to its termina- 
tion point, at C > 0.494. 
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FIG. 7. The same as in Fig. 5, but for 13 = 2. This branch 
is always unstable (as is shown by the middle plot demon- 
strating the instability growth rate vs. C) in its range of 
existence, < C < 0.219. Examples of the solution displayed 
in the upper part of the figure are given for C = 0.1 and 
C = 0.219. The lower part shows the new solution family 
found past the termination point of the unstable branch. Ex- 
ample of the new solutions are given for C — 0.22 (stable, the 
upper row) and C — 0.498 (just prior to the termination of 
the new family, the lower row). The instability of this branch 
sets in at C « 0.494, i.e., very close to the termination point. 

In the case of /3 = 2, direct simulations show that the 

instability of the anti-symmetric branch gives rise to re- 
arrangement of the solution into a very regular breather 
shown in Fig. 8 for C = 0.1 and k = 0.1. 

FIG. 8. The development of the instability accounted 

for by the imaginary oigenfrequency (with the growth rate 
~ 0.45) of the anti-symmetric branch, in the case of C = 0.1, 
K = 0.1. The top panels pertain to t = 4 (left), t = 100 
(middle) and t = 120 (right). The latter two have again been 
chosen close to the points where the oscillating amplitude of 
the resultant breather attains its maximum and minimum, re- 
spectively. The instability sets in around t 40; no external 
perturbation was added to the initial condition in this case. 
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D. Solution families which are asymmetric in the 
anti-continuum limit 



Additional branches of the solutions may start in the 
AC limit from asymmetric configurations, provided that 
A is still larger, namely for A > 2k. In particular, such an 
extra branch can be initiated by the following AC-limit 
solution excited at a single site n = no (here, (3 = 0), cf. 
Eqs. (17) and (18): 



A ± VA^ - 4k2 



>)• 



(19) 
(20) 



An example of this solution for the upper sign in Eq. (19) 
is shown, for A = 2, k = 0.5 and 6 = 1, in Fig. 9. Such 
asymmetric branches may be stable for sufficiently weak 
coupling (in this case, for C < 0.204). but they eventu- 
ally become unstable, and disappear soon thereafter (at 
C w 0.213, in this case). 
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is very close to the family's fixed propagation constant, 
A = 2. 




FIG. 10. The same as in Fig. 9, but generated by the 
expression (19) with the lower sign. This solution is always 
stable until it terminates at C ~ 0.458. Examples of the solu- 
tion for C = 0.2 and C = 0.458 (the latter case is chosen just 
prior to the termination of the branch) are shown, as usual, by 
means of their profiles and linear stability eigenfrequencies. 
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FIG. 9. The solution branch generated, in the AC limit, 
by the asymmetric expression (19)-(20) with the upper sign. 
The notation is the same as in Fig. 2. Two examples of the 
solution are shown for C = 0.1 and C = 0.213. The most un- 
stable eigenvalue is shown, vs. C, in the bottom panel. The 
instability sets in at C « 0.204, and the branch terminates at 
C « 0.213. 

The evolution of the instability (for C > 0.204) for this 
asymmetric branch is strongly reminiscent of that shown 
in Fig. 3, resulting in a persistent breathing state. 

The branch that commences from the AC expression 
(19) with the lower sign is shown for A = 2, k = 0.75 and 
^ = 1 in Fig. 10. The branch remains stable as long as 
it exists, i.e., for C < 0.46. At this point, it disappears 
colliding with the phonon band, whose upper edge is lo- 
cated, according to Eq. (10), at Aedge ~ 1.987, which 



III. GAP SOLITONS 
A. Solitons in the inner layer of the gap 

All the solutions that were examined in the previous 
section had their propagation constant above the upper 
edge of the phonon spectrum. Another issue of obvious 
interest is to study possible gap solitons (GSs), whose 
propagation constant is located inside the gap (11), i.e., 
below the lower edge of the phonon band. Unlike the soli- 
tons found above the band, GSs may persist up to the 
continuum limit. 

An example of such a solution for A = 0.75, k = —1, 
S = 0.9, and P — is shown in Fig. 11. In the AC limit, 
this branch starts with the expression (17). The branch 
is stable for small C, but then it becomes unstable due 
to oscillatory instabilities. The first two instabilities oc- 
cur at C = 0.242 and C = 0.349, as is shown in Fig. 
11. Past the onset of the instabilities, this branch con- 
tinues to exist (as an unstable one) indefinitely with the 
increase of C, and carries over into an (unstable) GS in 
the continuum limit. At large values of C, the distinct 
phonon bands, which are clearly seen in the example of 
the eigenvalue spectrum shown for C = 0.4 in Fig. 11, 
eventually collide and, due to their opposite Krein signs 
(see the definition and discussion of these in Ref. [16]), 
which gives rise to a whole set of oscillatory instabilities. 
The result is clearly seen in the example of the eigenvalue 
spectrum shown in the bottom panel of Fig. 11 for a large 
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value of the coupling constant, C = 4. The characteris- 
tic size of the instability growth rate (largest imaginary 
part of the eigenvalue) is nearly the same for C = 0.4 
and C = 4. in the latter case it being w 0.09. Notice, 
however, that, as the continuum limit is approached, the 
instabilities may be suppressed, in a part or completely, 
by finite-size effects (for an example of such finite-size 
restabilization, see Ref. [17]). 




S 0- 




FIG. 11. The branch of the gap-soHton solutions with 
A = 0.75, K = -1, 6 = 0.9, and ^ = 0. The upper part 

of the figure shows the norms of the two components of the 
soliton, and examples of the solutions for C = 0.1 (stable) 
and C = 0.4 (after the onset of the first oscillatory instabil- 
ity). The middle panel shows the instability growth rates, 
while the lower part of the figure gives an example of a solu- 
tion belonging to this branch, found at a much larger value of 
the coupling constant, C = 4. This solution family extends, 
as an unstable one, up to the continuum limit. 

The development of the oscillatory instability of GS 
belonging to the inner layer is displayed, for C > 0.242, 
in Fig. 12 for C = 0.4. 5 = 0.9, k = -1 and /? = 0. In 
this particular case, there are two oscillatory instabilities 
whose growth rates are in the interval 0.1 < < 0.2. As 
a result, symmetry breaking occurs, resulting in a shift of 
the central position of the soliton (from the site n = 50 
to n = 49). Oscillatory features in the dynamics are also 
observed in the latter case, and a small amount of energy 
is emitted as radiation. 
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FIG. 12. The evolution of the unstable gap soliton belong- 
ing to the inner layer, for C = 0.4, 5 = 0.9, k = — 1, and 
(3 = 0. The top panels shows the wave field distribution at 
t = 4 (left), t = 124 (middle) and t = 132 (right). Symme- 
try-breaking effects are clearly visible. A random perturba- 
tion of an amplitude 10^^ was added to the initial condition 
in order to catalyze the onset of the instability, which occurs 
at t > 40. 

Similar results were obtained for smaller values of A, 
for instance, A = 0.25. It was verified too that this sce- 
nario persists in the presence of the XPM nonlinearity 
(i.e., for P = 2), as it is shown in Fig. 13. In the latter 
case, the evolution of the instability with the increase of 
C is quite interesting, as it is nonmonotonic. The insta- 
bility first arises at C « 0.11 (due to a collision between 
discrete eigenvalues with opposite Krein signs). Subse- 
quent restabilization takes place at C « 0.16, but the 
solutions are unstable again for C > 0.36, and remain 
unstable thereafter, up to the continuum limit. 
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FIG. 13. The same as the previous figure, but for (5 = 2. 

Examples of the solutions arc shown for C = 0.2 (upper row, 
stable) and C = 1.6 (lower row, unstable due to several of 
oscillatory instabilities). The bottom panel demonstrates the 
nonmonotonic evolution of the instability of this solution with 
the increase of C. 

In the case of /3 = 2, the dynamical development of the 
oscillatory instabilities is similar to the (3 = case, again 
demonstrating symmetry-breaking effects. 



B. Solutions in the outer layer of the gap 

In all the cases considered in the previous subsections, 
the soliton's propagation constant A belonged to the in- 
ner layer of the gap, see Eq. (12). We have also examined 
the situation when A belongs to the outer layer defined 
in Eq. (12) (the outer layer exists unless S = 1). An 
example is shown in Fig. 14, where /3 = 0, k = — 1, 
S = 0.1, and A « 0.787 is chosen to be in the middle of 
the outer layer. In this case, we typically obtained de- 
localized solitons, sitting on top of a finite background 
(they are sometimes called "antidark" solitons). As can 
be observed from Fig. 14, such solutions may be sta- 



ble for sufficiently weak coupling, but become unstable 
as the continuum limit is approached, although they do 
not disappear in this limit (in Ref. [9] such delocalized 
solitons were found in the continuum counterpart of the 
present model). 
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FIG. 14. Solutions with non-decaying oscillatory back- 
ground for propagation constants belonging to the outer layer 
defined by Eq. (12). The top panel shows a stable solution 
for C = 0.158, and the bottom panel shows an unstable one 
of C = 0.576. These solutions are unstable for all C > 0.341. 



The instability development in the case of the outer- 
layer GSs is demonstrated, for C ~ 0.549, (3 ~ 0, k = —1 
and 6 = 0.1, in Fig. 15. In this case, the non-vanishing 
background is also perturbed by the instability, resulting 
in, plausibly, chaotic oscillations throughout the lattice. 
Symmetry-breaking effects, which shift the central peak 
from its original position, are observed too in this case. 
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FIG. 15. The time evolution of the outer-layer gap soli- 
tons for C = 0.549, /3 = 0, « = -1 and 5 = 0.1. The initial 

condition is perturbed by a random uniformly distributed per- 
turbation of an amplitude 10^^. The result of the instability 
is the excitation of background oscillations, as well as a shift 
of the soliton's peak from its original position. The top left 
and right spatial profiles correspond to t = 4 and t = 200, 
respectively. 



C. Stabilization of the gap solitons by mismatch 



The above considerations show that, inside the inner 
layer of the gap, it is easy to identify families of soliton 
solutions that persist in the continuum limit as C ^ oo. 
However, all the examples considered above showed that 
the solutions get destabilized at finite C and remain un- 
stable with the subsequent increase of C. Therefore, 
a challenging problem is to find solution families that 
would remain stable for large values of C. 

In fact, the introduction of a finite mismatch q (re- 
call it was set equal to zero in all the examples consid- 
ered above) may easily stabilize the discrete GSs. To this 
end, we pick up a typical example, with C = 0.5, k = — 1. 
6 = 0.5, A = 0.75, and /? = 0, when the GS exists but is 
definitely unstable in the absence of the mismatch. Fig- 
ures 16 and 18 show the effect of positive and negative 
values of the mismatch on the solitons. As is seen, large 
values of the positive mismatch can make the instability 
very weak, but cannot completely eliminate it. How- 
ever, sufficiently large negative mismatch readily makes 
the solitons truly stable. Thus, adding the negative mis- 
match is the simplest way to stabilize the solitons at large 
C, which is not surprising, as Eq. (11) demonstrates that 
the negative mismatch makes the gap broader. 





FIG. 16. A family of the gap-soliton solutions obtained for 

fixed values C = 0.5, k = -1, 5 = 0.5, A = 0.75, and = 0, 
by continuation to positive values of the mismatch parameter 
q. The top panel and the one beneath it show the evolution of 
the norms of the two components of the solution, and of the 
largest instability growth rate, with the increase of q. Other 
panels show examples of the solution (as usual, in terms of 
profiles of the two components and linear stability eigenval- 
ues) for g = 0, 9 = 2, and q = 5, from top to bottom. 



As an example of the dynamical evolution of unstable 
solitons in the case of positive mismatch, in Fig. 17 we 
display the case of C = 0.5, k = —1, /3 = 0, 5 = 0.5 
and q= I. In this case, the evolution leads to the estab- 
lishment of a breather with a rather complex dynamical 
behavior. 
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FIG. 17. The dynamical evolution in the unstable case with 
C = 0.5, K — —1, = 0, 5 = 0.5 and q = 1 (positive mis- 
match). The top left and right panels correspond to t = 4 
and t = 396, respectively. A complex pattern of the ampli- 
tude evolution is observed in this case. The initial condition 
contains a random perturbation with an amplitude 5 x 10~®. 
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FIG. 18. The same as in Fig. 16, but for negative val- 
ues of the mismatch. Examples of the solutions are given for 
(/ = 0, g = —2, and q = —5, from top to bottom. A random 
perturbation of an amplitude lO""* was used in this case. 

The instability development in the case of negative 
mismatch, q = —2.5, is demonstrated in Fig. 19. A 

localized breather with quasi-periodic intrinsic dynamics 
is observed in this case as an eventual state. 
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FIG. 19. The evolution of the unstable discrete gap soliton 
in with C = 0.5, k = -1, (3 = 0, 5 = 0.5 and q = -2.5 
(negative mismatch). The top left and right panels show the 
field profiles for t = 4 and t = 196, respectively. The time 
evolution of the amplitudes is shown in the bottom panel. 



D. Discrete counterparts of gap solitons from the 
Bragg-grating model 



As was shown in the introduction, the particular case 
of Eqs. (2) and (3) with k = 1, S = 1, and = 2 may be 
interpreted, with regard to the transformation (8), as a 
discretization of the standard Bragg-grating system (1). 
This continuum model gives rise to a family of exact GS 
solutions [19], 

^ = U{x) exp {-itcos6) , $ = V{x) exp {-it cos 9) , 

(21) 



U{x) = ^ sech (^x sine-^e^ ,V = -U*, (22) 

where the real parameter 6 takes values < ^ < tt. A 
part of this interval, < 6* < 6'cr ~ 1-01 (7r/2), is filled 
with stable solitons [18], while the remaining part con- 
tains unstable ones. 

All the discrete GSs considered above are not coun- 
terparts of the continuum solitons given by Eqs. (21) 
and (22). Establishing a direct correspondence between 
the latter ones and discrete solitons of Eqs. (2) and (3) 
is complicated by two problems: the transformation (8) 
does not have a continuum limit, and real symmetric or 
anti-symmetric GSs with |A| < k do not exist in the AC 
limit, as is seen from Eqs. (17) and (18), i.e., the usual 
starting point of the analysis is not available in this case. 

We have considered the discrete analogs of the Bragg- 
grating GSs in the following way. First, we took a formal 
discrete counterpart of the waveforms (21) and (22), and 
used them to obtain exact solutions of the (formal) dis- 
cretization of the Bragg-grating model of Eq. (1). Then 
the transformation (8) was applied to these solutions, 
and the thus obtained expressions were used as an initial 
guess for finding a numerically exact stationary solution 
of Eqs. (2)- (3). This procedure naturally generates new 
solitons, a crucial difference of which from all the types 
considered above is that they are truly complex solutions, 
see examples in Figs. 20 and 21. In this case, we have 
used 6 = -k/A. 
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FIG. 20. Absolute values of the fields in a new complex dis- 
crete gap-soliton solution, which was obtained, at C = 1.24, 
for the formal discretization of the Bragg-grating gap-soliton 
model. Then, the transformation (8) was applied to this solu- 
tion and it was used as an initial condition of the Newton for 
finding a solution of Eqs. (2)-(3). The continuous and dashed 
lines (which completely oevrlap) show the resulting profiles 
generated by the above-mentioned procedure, one which is 
an exact solution for the discretization of the Bragg-grating 
model and one which is an exact solution (identical in norm 
due to the nature of (8)) of Eqs. (2)-(3). 



FIG. 21. Left panels show profiles of real and imaginary 
parts of the fields ipn and </>„ in the discrete counterpart of 
the Bragg-grating gap soliton from Fig. 20. Right panels show 
stability eigenfrequencies for the same soliton. The upper and 
lower parts of the figure pertain to the soliton at C = 1.24 
(the same value as in Fig. 20), and to its continuation to 
the anti-continuum limit, C = 0. In the left panels, the cir- 
cles (joined by solid lines) refer to the real parts, while the 
stars (connected by dashed lines) correspond to the imaginary 
parts of the corresponding shown fields. 

Finally, linear stability eigenvalues were calculated for 
this new branch of the discrete GSs. The result (see Fig. 
22) is that this branch is unstable for all finite values of 
C, getting asymptotically stable in both limits C ^ 
and C CO (large values of C are not shown in Fig. 22); 
the stability regained in the latter limit complies with the 
above-mentioned finding that a subfamily of the contin- 
uum Bragg-grating gap solitons are dynamically stable. 
Notice that the natural norm of the contimmm soliton 
differs from that of the discrete one, given by Eq. (14), 
by an additional multiplier C^^/^ (which is proportional 
to the effective lattice spacing). Wc have chocked that 
the thus renormalized norm of the soliton converges as 
C ^ oo, although data for large C is not displayed in 
Fig. 22. 



Then, the solution was numerically continued, decreas- 
ing C, back to the AC limit, in order to identify its AC 
"stem". The result is shown in Fig. 21. Obviously, this 
AC state is very different from all those considered above 
(in particular, it is complex). 
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FIG. 22. The norm of the discrete counterpart of the 
Bragg-grating gap soliton of the model of Eqs. (2)- (3) (up- 
per panel), and its two unstable (imaginary) eigenfrequencies 
(lower panel shown by solid and dashed lines respectively), 
vs. the coupling constant. The continuation of the figure to 
larger values of C shows that the soliton becomes asymptoti- 
cally stable as C — > oo. 



IV. SOLITONS IN THE MODEL WITH THE 
QUADRATIC NONLINEARITY 

Stationary solutions of the SHG system (4) and (5) are 
looked for in an obvious form, cf. Eqs. (13): 



(23) 
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and in this case wc only consider the (most character- 
istic) case ^ = 1. The hnearization of Eqs. (4) and 
(5) demonstrates that one may expect termination of 
a sohton-sohition branch, due to its coUision with the 
phonon band, at (or close to) the point 



A = k/4 + C, 



and the gap between two phonon bands is 



(24) 



< A < k/4 



(25) 



(it exists only if k > 0). 

Stationary solutions were constructed, again, by means 
of continuation starting from the AC limit, where the ex- 
citation localized on a single site of the lattice assumes 
the form 



A, 



(26) 
(27) 



Note that solutions with the propagation constant be- 
longing to the gap (25) do not exist close to the AC limit. 
Indeed, the AC expression (27) shows that a necessary 
condition for its existence is 4A > k. On the other hand, 
A stays in the gap (25) if 4A < k, so the two conditions 
are incompatible. 

A typical example of a numerically found soliton 
branch is shown, for the value of the mismatch param- 
eter K = 0.9, in Fig. 23. This solution family is fixed 
by choosing A = 0.25, and starting from the expressions 
(26) and (27) with the upper sign. It is seen that the 
branch is stable for C < 0.015, but then it becomes un- 
stable, and eventually terminates at C w 0.025, in very 
good agreement with the prediction of Eq. (24). 



FIG. 23. The family of the discrete SHG solitons found for 
K = 0.9, 5=1, and A = 0.25. The top panel shows the 
evolution of the norms of the fundamental- (circles) and sec- 
ond-frequency (stars) components of the soliton, which are 
defined the same way as in Eq. (14), with the increase of C. 
Examples of solutions are displayed for C — 0.01, below the 
instability threshold, which is Ccr = 0.015 (the upper row), 
and for C = 0.024, just prior to the termination of the branch 
at C = 0.025 (the lower row). The bottom panel shows the 
evolution of the instability growth rate (imaginary part of the 
most unstable eigenvalue). 



Development of the instability of this SHG soliton 
branch for C > 0.015 was numerically examined through 
direct simulations, results of which are presented in Fig. 
24, for the case of C = 0.02, the corresponding instabil- 
ity growth rate being « 0.03. A random uniformly dis- 
tributed perturbation of an amplitude lO"** was added 
to the initial condition to accelerate the onset of the in- 
stability. The eventual result of the instability is the 
appearance of a breather with very regular periodic in- 
trinsic vibrations. 
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FIG. 24. Simulations of the instability development for the 
discrete SHG soliton in the case with C = 0.02, k = 0.9 and 
(5=1. The circles in the spatial profiles of the top panels 
denote the fundamental, and the stars denote the second har- 
monic. The left panel corresponds to t = 4, the middle to 
t = 660 (close to the point where the oscillating amplitude 
attains its maximum), and the right to 760 (close to a min- 
imum-amplitude point). The bottom panel shows the oscil- 
latory behavior at the central site for the fundamental (solid 
line) and second harmonics (dashed line). 



A similar result for the case without a gap in the 
phonon spectrum [see Eq. (25)] is shown in Fig. 25 
for K = —0.1, the other parameters being the same as in 
the previous case. This time, the branch terminates at 
C w 0.275, once again in complete agreement with the 
prediction of Eq. (24). 




FIG. 25. The same as in the upper part of Fig. 24, but for 
K = —0.1. The middle panel shows an example of the solu- 
tion for C = 0.1, and the bottom panel shows an example for 
C = 0.274, just prior to the termination of the branch (which 
happens at C = 0.275). 
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FIG. 26. The same as in Fig. 23, but for the lower sign in 
Eq. (27). The middle panel of the top subplot displays an 
example of a stable solution for C — 0.01, while the bottom 
panel shows a solution for C = 0.061, close to the termination 
of the branch. In this case, k = 0.75, 5=1, and A = 0.25. 



Lastly, another characteristic branch of solutions can 
be constructed starting from the pattern given by Eqs. 
(27) with the lower sign. This solution family is dis- 
played in Fig. 26, for k = 0.75, 5 = 1, and A = 0.25. 
The branch is stable for sufficiently weak coupling, but 
then it becomes unstable for C > 0.046. The branch dis- 
appears colliding with the phonon band at C « 0.062, 
once again in full agreement with the prediction of Eq. 
(24). 



An example of the development of instability of the 
present solution, that takes place at C > 0.046, is shown 
in Fig. 27 for C = 0.055 {k = 0.75; 6 = 1). A random 
initial perturbation with an amplitude 10~^ was added 
to the initial condition in this case. As is seen, the evo- 
lution results in complete destruction of the pulse into 
small-amplitude radiation waves. 
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FIG. 27. The evolution of the unstable SHG soliton in the 
case of C = 0.055, k = 0.75 and 6 = 1. The top left, middle, 

and right panels correspond to the profiles at t — 4, t = 160, 
and t = 400. The bottom panel, as before, shows the time 
evolution of amplitudes at the central site. 



V. CONCLUSION 

In this work, we have introduced a model which in- 
cludes two nonlinear dynamical chains with linear and 
nonlinear couplings between them, and opposite signs of 
the discrete diffraction inside the chains. In the case of 
the cubic nonlincarity, the model finds two distinct in- 
terpretations in terms of nonlinear optical waveguide ar- 
rays, based on the diffraction- management concept. A 
continuum limit of the model is tantamount to a dual- 
core nonlinear optical fiber with opposite signs of dis- 
persion in the two cores. Simultaneously, the system 
is equivalent to a formal discretization of the standard 
model of Bragg-grating solitons. A straightforward dis- 
crete second-harmonic-gencration [x*^^-*] model, with op- 
posite signs of the diffractions at the fundamental and 
second harmonics, was introduced too. Starting from 
the anti-continuum (AC) limit and gradually increasing 
the coupling constant, soliton solutions in the x^''' model 
were found, both above the phonon band and inside the 
gap. Above the gap, the solitons may be stable as long 
as they exist, but with transition to the continuum limit 
they inevitably disappear. On the contrary, solitons in 
the gap persist all the way up to the continuum limit. 
In the zero-mismatch case, they always become unstable 
before reaching the continuum limit, but finite mismatch 
may strongly stabilize them. A separate procedure had 
to be developed to search for discrete counterparts of the 
well-known Bragg-grating gap solitons. As a result, it 
was found that discrete solitons of this type exist at all 
values of the coupling constant C, but they appear to be 
stable solely in the limit cases C — and C = oo. Soli- 



In the cases when the solitons were found to be unsta- 
ble, simulations of their dynamical evolution reveal a va- 
riety of different scenarios. These include establishment 
of localized breathers featuring periodic, quasi-periodic, 
or very complex intrinsic dynamics, or destruction of one 
component of the soliton, as well as symmetry-breaking 
effects, and even complete decay of both components into 
small-amplitude radiation. The outcome depends on the 
type of the nonlinearity (cubic or quadratic), and on the 
nature of the unstable solution. 
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tons were also found in the x^^^ model. They too start 
as stable solutions, but then lose their stability. 
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